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1.0  Introduction 

Composites  of  fiber-reinforced  polymers  (FRPs)  respond  to  impact  loading  in  a  variety  of 
inelastic  modes  of  deformation,  depending  on  such  variables  as  the  thickness  of  the  composite 
section  relative  to  the  width  of  the  impactor,  the  mechanical  properties  of  the  fibers  and  the 
matrix  material,  the  bonding  between  the  fibers  and  the  matrix,  the  weave  of  the  fibers,  and  the 
impact  velocity.  In  an  effort  to  understand  these  inelastic  modes,  much  research  has  recently 
focused  on  the  mechanics  of  composites  under  impact  at  the  meso-scale.  This  research  includes 
some  numerical  simulations  of  composites  at  the  meso-scale,  in  which  the  matrix  material  and 
yarns  of  fibers  are  modeled  by  separate  finite  elements.  Such  simulations  shed  light  on  the 
inelastic  modes  of  interest,  but  they  generally  require  too  much  mesh  refinement  to  be  useful  at 
the  macroscopic  scale  of  design  analysis.  Therefore,  in  the  analysis  of  FRP  applications,  some 
degree  of  material  homogenization  is  generally  necessary. 

Two  approaches  are  now  available  to  model  FRP  composites  in  EPIC  [1].  In  one  approach, 
homogenization  is  restricted  to  the  individual  materials.  With  this  restriction,  a  uniaxial-stress 
model  for  the  fibers  is  implemented  in  bar  elements,  while  an  isotropic  elastic-plastic  model  for 
the  matrix  is  implemented  in  3D  solid  elements,  and  the  two  element  types  interact  through 
common  nodes  in  the  mesh  [2],  A  sketch  of  this  modeling  approach  is  shown  in  Figure  1. 
Homogenization  refers  to  the  representation  of  multiple  yams,  including  those  from  multiple 
plies  of  fabric,  by  the  same  bar  element.  Similarly,  the  matrix  material  separated  by  plies  of 
fabric  is  represented  by  the  same  solid  element.  The  primary  advantage  of  this  approach  is  that 
the  mechanical  properties  of  the  fabric  can  be  acquired  from  laboratory  experiments  and  applied 
directly  to  the  bar  elements.  The  primary  disadvantage  is  that  the  bar  elements  occupy  no 
volume,  the  solid  elements  necessarily  occupy  one-hundred  percent  of  the  composite  volume, 
and  the  properties  of  the  solid  elements  therefore  cannot  be  acquired  from  experiments;  they 
must  be  extrapolated  from  the  properties  of  the  matrix  material. 

Composite  layup  Composite  FE  model 


Figure  1.  One  approach  to  modeling  composites  in  EPIC:  bar  elements 
represent  the  yarns  of  fibers  and  solid  elements  represent  the  matrix. 

In  the  second  approach  to  modeling  composites  in  EPIC,  the  responses  of  the  fibers  and  the 
matrix  material  are  homogenized  into  a  single  material  model,  and  that  model  is  implemented  in 
3D  solid  elements.  This  approach  represents  full  homogenization  of  the  composite  material,  and 
it  has  been  implemented  via  various  anisotropic  material  models  [3-6]  in  other  finite-element 
codes.  The  primary  advantage  of  this  approach  is  that  the  difference  between  modeling  an 
isotropic  metal  and  an  orthotropic  composite  is  confined  to  the  material  model,  so  that  mesh 
generation  is  identical  for  all  materials.  The  primary  disadvantage  is  the  relative  complexity  of 
the  material  model  that  is  required  to  adequately  represent  an  orthotropic  continuum.  This  report 
documents  the  integration  of  orthotropic  plasticity  and  failure  models  into  the  existing 
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orthotropic  elasticity  model  in  EPIC  for  the  purpose  of  modeling  FRP  composites  under  ballistic 
impact  with  full  material  homogenization.  It  includes  a  section  of  example  computations  to 
demonstrate  the  orthotropic  models  and  evaluate  how  well  they  apply  to  FRP  composites. 
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2.0  Orthotropic  Model  Formulation 

Several  models  have  been  proposed  in  the  literature  for  specific  fiber-reinforced  polymers. 
In  general,  these  models  include  assumptions  germane  to  the  specific  composites  of  interest,  and 
generate  computed  results  that  correlate  well  with  the  laboratory  tests  under  consideration.  For 
example,  visco-plastic  models  for  glass-fiber  composites  are  typically  simplified  by  only 
admitting  elastic  strains  in  the  directions  of  the  fibers  [5].  This  simplifying  assumption  is 
supported  by  the  linearity  of  the  stress-strain  curve  up  to  failure  in  tensile  tests  along  the 
directions  of  the  fibers,  especially  in  comparison  to  the  non-linear  stress-strain  curves  obtained 
from  off-axis  tensile  tests. 

A  somewhat  more  general  approach  has  been  taken  for  EPIC,  based  on  the  desire  to  treat  as 
many  composite  materials  as  practical  with  one  model,  but  also  based  on  the  recognition  that  the 
complexities  of  composite  materials  make  an  optimal  model  difficult  to  identify  before 
exhaustive  tests  and  computations  have  been  perfonned.  The  approach  consists  of  adding 
ortho  tropic  plasticity  and  failure  models  to  the  existing  ortho  tropic  elasticity  model  in  EPIC. 
This  section  describes  the  implementation  of  those  models. 

2.1  Material  Reference  Frame 

When  finite  rotations  or  large  displacements  are  considered,  a  finite-element  analysis  must 
account  for  the  change  in  orientation  of  the  material  with  time.  If  the  material  is  isotropic,  the 
constitutive  models  are  often  cast  in  rate  fonn,  with  a  suitable  stress  rate  to  update  the  Cauchy 
stress  tensor  in  the  rotating  material.  But  when  the  material  is  anisotropic,  it  also  becomes 
necessary  to  monitor  the  orientation  of  the  principal  material  directions,  since  they  may  each 
exhibit  unique  behavior.  This  can  be  done  most  naturally  by  evaluating  the  constitutive 
equations  in  a  frame  that  rotates  with  the  material,  i.e.,  a  material  reference  frame. 

The  equations  of  motion  must  be  updated  in  a  global,  inertial  reference  frame.  If  the 
constitutive  equations  are  evaluated  in  the  material  frame,  it  becomes  necessary  to  transfonn 
between  the  two  frames  each  timestep.  Figure  2  depicts  the  two  reference  frames,  and  the 
transformations  necessary  to  evaluate  the  constitutive  models  in  a  material  reference  frame.  In 
the  figure,  the  global  frame  is  represented  by  X,  Y  and  Z;  and  the  material  frame  is  represented 
by  X1  ,  X2  and  X3.  The  rotation  tensor,  R,  transforms  vectors  and  tensors  from  the  material 
frame  to  the  global  frame.  It  is  composed  of  the  initial  orientation  of  the  global  frame  with 
respect  to  the  material  frame,  Rt=o >  and  the  subsequent  rotation  of  the  material  during  the 
computation,  Rt>0.  The  algorithm  in  EPIC  for  updating  the  rotation  tensor  was  first  proposed  by 
Dienes  [7].  Rather  than  referring  to  the  initial  configuration  each  timestep,  the  algorithm  updates 
the  rotations  using  the  rates  of  deformation.  Anderson  et  al  [8]  discuss  the  algorithm  in  detail. 

The  practice  of  updating  the  Cauchy  stress  tensor  of  an  orthotropic  material  in  the  material 
frame  exactly  accounts  for  finite  rotations,  but  does  not  account  for  finite  shearing  defonnations. 
This  practice  will  provide  adequate  accuracy  when  shear  strains  are  small,  with  errors  accruing 
only  when  non-negligible  material  strength  exists  under  finite  shear  strains.  A  total  Lagrangian 
formulation  would  account  for  finite  shear  strains  because  the  stress  update  is  perfonned  on  the 
Piola-Kirchoff  stress  tensors,  which  are  cast  with  respect  to  the  initial  unstrained  configuration. 
However,  the  advantages  of  such  a  formulation  are  not  clear,  given  the  current  limitations  in 
accuracy  of  homogenized  inelastic  models  of  composites. 
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Material  orientation  ^  =  Rt= 0+  ^t> 0  Current  orientation  (t  >  0) 

-  define  material  parameters  -  update  transient  rotations,  Rt>0 

-  evaluate  constitutive  models  -  evaluate  equations  of  motion 

(update  stresses) 

Figure  2.  The  rotations  that  transform  an  element  of  material  between 
its  initial,  material,  and  current  orientations. 

2.2  Orthotropic  Elasticity 

The  orthotropic  elasticity  model  in  EPIC  [1]  is  a  generalization  of  small-strain  orthotropic 
elasticity  to  finite  strains  by  treating  each  timestep  At  of  the  explicit  computation  as  an  increment 
of  small  strain.  As  a  result,  the  incremental  stress-strain  relationship  remains  linear,  and  takes 
the  following  form. 

De=—C1:  A<t  (1) 

At  v  ' 

In  this  equation,  De  is  the  elastic  rate-of-deformation  tensor  that  results  from  an  additive 
decomposition  of  the  rate-of-defonnation  tensor  into  elastic  and  plastic  parts,  D  =  De  +  Dp  ;  Act 
is  the  increment  in  Cauchy  stress;  and  C  is  the  tensor  of  elastic  moduli.  In  EPIC,  this  equation  is 
implemented  via  central  differences  because  the  velocity  gradients  of  the  rate-of-defonnation 
tensor  are  evaluated  at  the  midsteps.  Expressing  this  incremental  elastic  relationship  in  matrix 
form, 
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In  this  expression,  the  three  Young’s  moduli  ( E l5  E2,  E 3)  and  the  six  Poisson’s  ratios  (v12, 
v2i,  iq3,  V31,  v23,  v32)  can  be  obtained  from  tensile  tests  in  the  three  principal  material 
directions.  But  they  are  not  all  independent;  symmetry  considerations  impose  three  conditions 
on  the  elastic  moduli  that  reduce  the  number  of  independent  constants  from  twelve  to  nine: 
v2i  —  v12E2/E1,  v3i  =  V13E3/E1  and  v32  =  v23£,3/£'2. 

Because  the  stresses  are  updated  each  timestep  from  the  deformations,  the  incremental 
elastic  relations  must  be  inverted, 

A  a  =  At  C.  De  (3) 

Or,  in  matrix  notation, 
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where /?  — (1  v12v21  v23v32  v31v13  2v21v32v13)/(£'1£'2£'3). 


2.3  Orthotropic  Plasticity 


The  primary  limitation  on  the  generality  of  the  plasticity  model  is  the  choice  of  the 
functional  form  of  the  yield  surface.  The  Hill  yield  surface  [9]  was  originally  introduced  to 
model  orthotropic  metal  plasticity.  Several  other  orthotropic  yield  surfaces  have  been  proposed 
since  then,  with  the  aim  to  improve  accuracy  for  specific  material  classes,  but  no  consensus  has 
been  reached  on  the  best  functional  fonn  for  FRP  composites.  As  a  result,  the  Hill  yield  surface 
has  been  implemented  in  EPIC,  with  the  goal  of  fully  evaluating  its  applicability  to  FRP 
composites  as  test  data  become  available.  This  yield  surface  is  quadratic  in  the  stress 
components, 

1 

/ O)  =  2  \E&22  -  03b)2  +  GO33  -  <Al)2  +  #Oll  -  °22)2]  +  NOi22  +  Mo ?3  +  Lo f3  -  O2  =  0 

(5) 


In  this  equation,  the  subscripts  on  the  stresses  refer  to  the  orthogonal  principal  directions  of 
the  material.  For  a  fiber  composite  composed  of  0/90-degree  plies,  two  of  the  principal 
directions  align  with  the  fiber  directions,  and  the  third  principal  direction  is  nonnal  to  the  plane 
of  the  composite. 


The  model  parameters  are  F,  G,  H,  N,  M  and  L;  and  the  average  radius  of  the  yield  surface  in 
stress  space  is  represented  by  a.  The  normal  yield  stresses  (cr^,  (iyZ2,  and  ay3)  and  shear  yield 
stresses  (ay2,  ay3,  and  o23)  are  aligned  with  the  principal  material  directions,  and  they  are  used  to 
determine  the  model  parameters: 


F  = 


(6) 
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G  = 


H  = 


N  = 


M  = 


L  = 


(7) 

(8) 

(9) 

(10) 

(11) 


The  yield  stresses  are  the  inputs  to  EPIC  that  define  the  Hill  yield  surface,  and  the  model 
parameters  are  computed  from  the  yield  stresses  according  to  equations  6-11. 

Material  hardening  with  inelastic  strains  has  been  implemented  through  the  average  normal 
radius  of  the  yield  surface,  o. 

<*  =  \  (°iyi  +  a22  +  ^33 )  ( 1  +  cif?)  (12) 

In  this  equation,  the  initial  average  yield  stress  is  factored  by  the  hardening  tenn,  (l  -I-  Cjfp2), 
£p  is  the  effective  plastic  strain,  and  C1  and  C2  are  the  hardening  parameters.  Since  £p  is  a 
scalar,  hardening  is  isotropic,  i.e.,  it  does  not  depend  upon  the  orientation  of  the  inelastic  strains. 
Hardening  in  this  form  can  be  interpreted  as  a  uniform  scaling  of  the  yield  surface  in  all 
directions  of  stress  space. 

To  complete  a  model  of  inelastic  defonnations,  the  inelastic  strains  must  be  characterized  by 
a  flow  rule.  An  associated  flow  rule  enforces  the  principal  of  normality,  associating  the  direction 
of  inelastic  rates  of  deformation  with  the  nonnal  to  the  yield  surface.  While  some  frictional 
materials  do  not  obey  the  principal  of  normality,  most  others  do.  As  a  result,  an  associated  flow 
rule  has  been  chosen  for  the  ortho  tropic  inelastic  model  in  EPIC, 

DP=A^  (13) 

da 

In  this  equation,  Dp  is  the  plastic  rate-of-defonnation  tensor;  is  the  nonnal  to  the  yield  surface 

/,  and  the  scalar  A  is  the  plastic  proportionality  factor.  When  the  yield  surface  is  quadratic  in  the 
stress  components,  like  the  Hill  yield  surface,  the  normal  to  the  surface  is  linear  in  the  stress 
components,  and  the  associated  flow  rule  takes  the  following  form, 

Dp  =  A  L:a  (14) 

In  this  equation,  L  is  a  tensor  of  constant  coefficients.  In  matrix  notation, 
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2.4  Elastic-Plastic  Constitutive  Update 

The  incremental  elastic-plastic  constitutive  update  is  expressed  by  substituting  De  —  D  — 
Dp  into  the  incremental  elastic  update  of  equation  3, 


Act  =  At  C:  (D  —  Dp)  (16) 

Further  substitution  of  the  flow  rule  in  equation  14  for  Dp  yields, 

Ao-  =  At  C:  (D  —  A  L:  a)  (17) 


The  stress  increment  is  then  expressed  explicitly  as  A<r  =  an+1  —  <jn,  in  which  an  and 
crn+1are  the  stress  tensors  at  the  beginning  and  end  of  the  (n+ 1  )th  timestep,  respectively.  Also, 
the  variable  is  substituted  for  the  product  A  At. 

an+1  —  an  +  At  C:  D  —  <fC:  L:  a  (18) 

Finally,  rearranging  terms  gives, 

an  +  AtC:D  —  an+1  +  %C\  L:  a  (19) 

The  two  terms  on  the  left  side  of  equation  19  are  known  at  the  beginning  of  the  elastic- 
plastic  constitutive  update,  and  are  referred  to  collectively  as  the  elastic  predictor  of  the  updated 
stress,  a*  =  an  +  At  C:  D.  The  first  term  on  the  right  side  is  the  desired  solution  for  the  updated 
stress,  and  it  is  limited  by  the  yield  surface.  The  second  tenn  on  the  right  side,  <fC:  L:  a,  is  the 
plastic  corrector,  since  it  represents  the  difference  between  the  elastic  predictor  and  the  updated 
stress.  The  plastic  corrector  is  linear  in  the  stress,  with  constant  coefficients  C  and  L,  and  the 
unknown  scalar  coefficient 

If  the  elastic  predictor  lies  within  the  yield  surface,  (<rn  +  At  C:  D)  —  a2  <  0  ,  then  the  rate 
of  deformation  is  entirely  elastic  for  that  timestep,  the  plastic  corrector  is  zero,  and  the  updated 
stress  is  set  equal  to  the  elastic  predictor.  If  the  elastic  predictor  exceeds  the  yield  surface,  then 
the  plastic  corrector  is  nonzero  and  the  expression  for  an+1  must  be  solved  simultaneously  with 

/On+i)  -  er2  =  0. 

The  method  of  this  solution  depends  upon  the  choice  of  a  in  the  plastic  corrector.  To 
maintain  consistency  with  the  governing  differential  equation,  a  must  be  on  the  interval  from 
an  to  <Tn+1 .  Two  schemes  of  first-order  accuracy  are  obtained  if  a  is  assigned  its  values  at  the 
endpoints  of  the  timestep.  If  the  value  is  at  the  beginning  of  the  timestep,  a  =  an,  then  the 
nonnal  to  the  yield  surface  during  the  timestep  is  approximated  by  its  value  at  the  beginning  of 

the  timestep,  i.e.,  an,  the  plastic  corrector  becomes  <fC:  L\  an,  and  the  resulting  finite- 

difference  scheme  is  referred  to  as  a  forward-Euler  scheme.  If,  on  the  other  hand,  the  value  is  at 
the  end  of  the  timestep,  a  —  <rn+1,  then  the  nonnal  to  the  yield  surface  is  approximated  by  its 
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value  at  the  end  of  the  timestep,  ^*11:  crn+1,  the  plastic  corrector  becomes  <fC:  L\  an+1,  and 

the  resulting  scheme  is  referred  to  as  a  backward-Euler  scheme.  In  developing  an  option  for 
orthotropic  elasto-plasticity  in  EPIC,  both  of  these  schemes  have  been  implemented.  The 
remainder  of  this  section  describes  the  implementations  of  these  schemes,  while  the  following 
section  on  numerical  examples  includes  a  comparison  of  their  performances. 

Figure  3  depicts  the  solution  procedure  in  stress  space  for  both  finite-difference  schemes. 
(Due  to  graphical  limitations,  the  six-dimensional  stress  space  is  reduced  to  two  dimensions.) 
The  yield  surface  is  represented  by  the  red  line,  and  denoted  by  /(< r)  —  a2  —  0.  The  state  of 
stress  at  the  beginning  of  the  timestep,  an  ,  is  located  inside  the  yield  surface,  while  the  elastic 
predictor,  a* ,  is  outside  the  yield  surface.  This  represents  the  most  general  case  of  an  elasto- 
plastic  stress  update. 


forward-Euler  scheme 


backward-Euler  scheme 


f(<r)  -  a2  =  0 


Figure  3.  Depictions  of  the  forward-Euler  (left)  and  backward-Euler 
(right)  finite-difference  schemes  in  time  for  the  elastic-plastic 
stress  update.  The  six-dimensional  stress  space  is  shown  in  two 
dimensions,  and  the  yield  surface  is  indicated  by  the  red  lines. 

On  the  left  side  of  the  figure  is  the  forward-Euler  scheme.  In  this  scheme,  it  is  first 
necessary  to  find  the  stress  state  along  the  path  of  the  elastic  predictor  that  lies  on  the  yield 
surface.  This  stress  state  is  denoted  by  <7^,  and  it  is  interpolated  between  the  known  stress  states 
an  and  a*  as  <7^  =  (1  —  a)  an  +  at 7*.  It  is  determined  by  solving  for  a  in  the  yield  condition, 


f{?yn)  -o2  =  /(( 1  -  a)  an  +  aa *)  -  cr2  =  0 


(20) 


Since  the  yield  surface  is  quadratic  in  the  stresses,  which  are  in  turn  linear  in  the  unknown 
a,  this  expression  reduces  to  a  quadratic  equation  in  a,  and  can  therefore  be  solved  explicitly. 
The  direction  of  the  plastic  corrector  is  then  calculated  from  <7^  as  L:  <7^,  and  it  is  represented 
in  the  figure  by  the  solid  arrow  emanating  from  <7^ .  The  expression  for  the  updated  stress  is  then 
substituted  into  the  equation  of  the  yield  surface, 


/On+i)  -  cr2  =  f(tF*  -  L.  a])  -  a2  =  0 


(21) 


And  this  equation  is  then  solved  for  the  unknown  scaling  factor,  which  in  turn  provides  the 
updated  stress  state,  <7n+1.  The  expression  for  also  reduces  to  a  quadratic  equation,  with  an 


8 


UNCLASSIFIED 


UNCLASSIFIED 


explicit  solution.  Therefore,  the  forward-Euler  scheme  for  the  elasto-plastic  stress  update  is 
explicit,  and  its  graphical  representation  is  the  connection  of  stress  states  a*  and  crn+1  by  the 
dashed  arrow  on  the  left  side  of  Fig  3. 

The  right  side  of  Fig.  3  depicts  the  backward-Euler  scheme.  In  this  scheme,  the  direction  of 
the  plastic  corrector,  (C:  L:  an+1 ,  depends  on  the  unknown  updated  stress,  crn+1,  which  lies  on 
the  yield  surface.  There  is  therefore  no  need  to  find  the  stress  state  along  the  path  of  the  elastic 
predictor  that  lies  on  the  yield  surface.  Substituting  <rn+1for  a  in  the  elasto-plastic  stress  update 
yields, 

on  +  AtC:D  =  (I  +  %C:L):an+1  (22) 


where  /  is  the  fourth-order  identity  tensor.  Substituting  a*  for  the  elastic  predictor  on  the  left 
side  of  equation  22,  and  solving  for  the  updated  stress, 

an+i  —  (/  +  ‘fC:  I*)- 1 :  a*  (23) 


Due  to  the  symmetry  of  the  stress  tensor,  this  expression  represents  a  system  of  six 
equations.  They  express  the  unknown  stresses  on+1  as  rational  algebraic  functions  of  the 
unknown  scalar  and  they  must  be  solved  simultaneously  with  the  yield  condition  f{an+1)  — 
a2  —  0.  Analytic  expressions  for  the  inverse  of  (/  +  (T:  L)  have  been  derived  for  EPIC,  and  a 
Newton-Raphson  iterative  method  on  the  yield  condition  has  been  implemented  for  the  solution 
of  the  updated  stresses. 


The  Newton-Raphson  method  begins  with  an  estimate  of  An  iterative  loop  is  then  entered 
in  which  an  estimate  of  the  updated  stress  is  computed  from  using  the  inverse  of  (/  +  <fC:  L). 
On  the  right  side  of  Fig.  3,  this  estimated  stress  state  is  denoted  by  <rln+1,  where  the  superscript  i 
indicates  the  iteration  count.  Since  this  stress  state  has  been  computed  from  an  estimate  of 
(Tln+1  does  not  lie  on  the  yield  surface  in  the  figure.  The  yield  condition,  (< Tln+1 )  —  a2  ,  and  its 


derivative  with  respect  to  <f, 


d/Ki+i)  _  5/  da 


are  then  computed  and  the  estimate  of  is 


d%  da  df 

improved  via  the  Newton-Raphson  method.  The  approximate  solution  to  the  updated  stress, 
°7i+i’  approaches  the  yield  surface  with  each  successive  iteration  on  the  value  of  f((Tln+i)  ~  °2 > 
until  the  desired  tolerance  on  the  root  of  this  expression  has  been  met.  The  final  value  of  oln+1  is 
then  assigned  to  <xn+1,  and  the  algorithm  exits  the  iterative  loop. 


The  quality  of  the  initial  estimate  of  affects  both  the  reliability  and  efficiency  of  the 
iterative  method.  The  reliability  is  reflected  by  the  rate  of  success  in  finding  the  solution,  and  the 
efficiency  is  reflected  by  the  number  of  iterations  required.  An  initial  estimate  of  that  has 
proven  very  reliable  and  efficient  is  obtained  by  scaling  the  elastic  predictor  back  to  the  yield 
surface  along  the  line  to  the  origin  of  stress  space,  computing  the  direction  of  the  plastic 
corrector  at  that  stress  state,  and  projecting  a*  —  an  onto  that  direction. 

When  plasticity  is  coupled  with  orthotropic  elasticity,  the  tensor  of  elastic  constants 
generally  rotates  the  plastic  corrector  away  from  the  nonnal  to  the  yield  surface,  AC:  as 

depicted  for  both  schemes  in  Fig.  3.  In  the  special  case  of  isotropic  elasticity,  the  elasto-plastic 
stress  update  can  be  expressed  in  terms  of  the  deviatoric  stress,  the  tensor  of  elastic  constants 
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reduces  to  a  scalar  multiple  of  the  identity  tensor,  C  —  2GI,  and  the  plastic  corrector  is  therefore 
normal  to  the  yield  surface. 

Although  the  forward  and  backward  Euler  schemes  are  both  first-order  accurate,  a 
significant  practical  difference  exists  between  the  two  methods,  and  it  results  from  representing 
the  normal  to  the  yield  surface  by  its  value  either  at  the  beginning  or  end  of  the  timestep.  If  the 
normal  to  the  yield  surface  is  represented  by  its  value  at  the  end  of  the  timestep,  as  in  the 
backward  Euler  scheme,  then  the  direction  of  the  plastic  corrector  varies  with  the  solution,  crn+1. 
Given  a  smooth  yield  surface,  this  variation  in  the  directions  of  the  plastic  correctors  allows  them 
to  span  the  entire  space  of  elastic  predictors,  guaranteeing  a  solution  to  the  backward  Euler 
scheme. 

If  the  normal  to  the  yield  surface  is  represented  by  its  value  at  the  beginning  of  the  timestep, 
as  in  the  forward-Euler  scheme,  then  the  direction  of  the  plastic  corrector  is  fixed  for  all  points 
on  the  yield  surface,  and  the  plastic  correctors  span  only  a  subset  of  the  space  of  elastic 
predictors.  An  elastic  predictor  that  lies  outside  the  space  of  forward-Euler  plastic  correctors  is 
depicted  in  Fig.  4.  Here,  the  fixed  direction  of  the  plastic  corrector  has  been  copied  to  several 
locations  along  the  yield  surface,  as  indicated  by  the  dashed  arrows,  and  the  resulting  space 
spanned  by  the  plastic  correctors  is  shaded  gray.  This  figure  demonstrates  the  loss  of  a  solution 
to  the  forward-Euler  scheme  that  may  occur  when  the  timestep  is  too  large  for  the  rate  of 
loading. 


Figure  4.  An  example  of  an  elastic  predictor  that  lies  outside  the  space 
(shaded  gray)  spanned  by  the  plastic  correctors  in  a  forward-Euler 
scheme.  In  this  scenario,  the  timestep  is  too  large  for  the  rate  of  loading. 

The  Hill  yield  function  reduces  to  the  von  Mises  yield  function  when  F  —  G  —  H  —  1  and 
N  —  M  —  L  —  3.  When  the  von  Mises  yield  function  is  implemented  with  isotropic  elasticity, 
the  initial  estimate  of  in  the  backward-Euler  scheme  provides  the  converged  solution,  and  the 
scheme  reduces  to  the  familiar  radial-return  algorithm. 

2.5  Orthotropic  Damage  and  Failure 

To  complete  an  orthotropic  material  model  for  the  ballistic  impact  of  composites,  an 
orthotropic  failure  model  has  also  been  implemented  in  EPIC.  This  model  acts  in  conjunction 
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with  the  orthotropic  elasticity  and  plasticity  models  described  in  the  preceding  sections.  It 
compiles  components  of  the  plastic  natural-strain  tensor  associated  with  the  principal  material 
directions,  and  fails  the  material  by  eliminating  the  deviatoric  stresses  when  any  of  the  plastic 
strain  components  reaches  its  user-supplied  critical  value. 

The  plastic  natural-strain  tensor  can  be  derived  from  the  rate-of-deformation  tensor, 
beginning  with  the  relation  (Malvern,  [11])  between  the  rate-of-deformation  tensor  D  and  the 
infinitesimal  line  segment  dx, 

Yt(ds2)  =  2dx-D-dx  (24) 

In  this  expression,  ds  is  the  magnitude  of  the  line  segment,  dx  —  ds  n,  where  n  is  a  unit 
vector  in  the  direction  of  dx.  This  expression  indicates  that  the  rate-of-defonnation  tensor 
measures  the  instantaneous  rate  of  change  of  the  squared  length  of  an  infinitesimal  line  segment. 
Since  ds  is  associated  with  the  spatial  segment  dx,  not  the  material  segment  dX,  integration  of 
this  expression  in  time  with  respect  to  an  inertial  reference  frame  does  not  generally  produce  a 
material  metric.  In  this  implementation,  however,  the  reference  frame  is  the  material  frame, 
which  moves  and  rotates  with  the  material,  and  temporal  integration  is  meaningful  in  the  small- 
strain  limit.  Therefore,  replacing  the  differentials  by  finite  differences  between  tn  and  tn+1, 

( dsn+1 )2  —  (dsn)2  =  2  A t  dx  -  D  ■  dx  (25) 

And  factoring  the  left  side, 

( dsn+1  +  dsn)(dsn+1  —  dsn )  =  2  A t  dx  -  D  ■  dx  (26) 


Denoting  the  mid-step  segment  magnitude  as  dsn+1/2  —  (dsn+1  +  dsn)/2  and  dividing  by  twice 
its  square, 

(dsn+1  -  dsn)/dsn+1/2  =  At  •  D  ■  (27) 

asn+ 1/2  asn+ 1/2 

For  the  increment  in  time  from  tn  to  tn+1,  the  two  fractions  on  the  right  side  are  the  unit  vector, 
n,  and  the  left  side  is  an  increment  in  natural  strain  As  along  dx. 

As  =  At  n  D  ■  n  (28) 


From  the  previous  section,  the  rate-of-defonnation  tensor  is  composed  of  elastic  and  plastic 
parts,  D  —  De  +  Dp .  Similarly,  the  natural-strain  increments  can  be  decomposed  into  elastic  and 
plastic  parts.  Substituting  these  decompositions  into  equation  28, 

Ase  +  Asp  =  At  n  ■  (De  +  Dp)  n  (29) 

When  the  deformation  is  entirely  elastic,  equation  29  must  reduce  to  A se  =  At  n  ■  De  ■  n. 
Subtraction  of  this  expression  for  A se  from  equation  29  gives, 

Aep  =  A  tn  Dp  n  (30) 

This  expression  defines  the  plastic  natural-strain  tensor  in  material  coordinates. 
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3.0  Numerical  Examples 


3.1  Taylor  Anvil 

Computations  of  a  cylinder  impacting  a  rigid  surface  demonstrate  the  effects  of  the 
ortho  tropic-plasticity  model  parameters.  They  also  provide  some  verification  of  the  model’s 
implementation  in  EPIC. 

The  top  left  comer  of  Fig.  5  shows  an  overhead  view  of  the  cylinder’s  mesh  in  the  initial 
configuration,  with  the  principal  material  directions  superimposed.  The  mesh  is  composed  of 
40,320  tetrahedral  elements.  The  cylinder  is  3.0  cm  in  length  and  1.5  cm  in  diameter,  and  its 
impact  velocity  is  150  m/s.  The  top  right  comer  of  the  figure  shows  the  results  of  the 
computation  using  the  isotropic  model  parameters  listed  in  Table  1.  The  results  are  colored  by 
effective  plastic  strain  at  100  ps  after  impact,  when  plastic  deformation  is  complete.  The  results 
of  the  isotropic  computation  are  identical  to  those  obtained  by  using  EPIC’s  model  for  metal 
plasticity,  which  employs  the  von  Mises  yield  surface  and  the  Johnson-Cook  strength  model.  The 
isotropic  computation  therefore  verifies  that  the  orthotropic  model  is  correctly  implemented  in 
the  special  case  of  material  isotropy. 

Table  1.  Isotropic  Material  Parameters  Used  for  the  Computation 
of  a  Cylinder  Impacting  a  Rigid  Surface 

p  —  7900  kg/m3 
E-i  —  E2  —  £3  —  192  GPa 
On  —  t/13  —  G23  —  80  GPa 

v12  —  v13  —  v23  —  0-2 
°ii  =  °22  =  °33  =  125  MPa 
ai2  =  ai3  ~  °23  =  101  MPa 

C1  =  0,  c2  =  1 

ppf  _  ppf  _  ppf  _  ppf  _  ppf  _  ppf  _  1 
_ fc11  fc33  fc12  fc13  t?.3  L _ 


By  comparison  with  the  orthotropic  computations,  the  isotropic  computation  also  helps  to 
detennine  the  effects  of  ortho tropy.  The  remaining  six  plots  in  Fig.  5  show  the  results  of  the 
orthotropic  computations,  each  with  the  same  material  parameters  as  the  isotropic  computation, 
except  a  50  percent  increase  in  one  of  the  yield  stresses:  ,  a^2,  er^3  ,  crf2>  °\3  or  °23-  Parts  c,  d 

and  e  show  the  plastic  strains  due  to  increases  in  the  nonnal  yield  strengths,  and  parts  f,  g  and  h 
show  the  plastic  strains  due  to  increases  in  the  shear  yield  strengths.  When  the  in-plane  nonnal 
strengths  are  increased  (parts  c  and  d),  or  the  out-of-plane  shear  strengths  are  increased  (parts  g 
and  h),  the  deformation  of  the  cylinder  is  reduced  in  the  respective  direction. 

These  computations  were  perfonned  with  both  the  forward-Euler  and  backward-Euler 
schemes  described  in  the  preceding  section.  To  maintain  stability,  the  forward-Euler  scheme 
required  a  timestep  about  one  tenth  that  of  the  Courant  condition  in  the  initial  configuration,  and 
about  one  third  that  of  the  Courant  condition  in  the  final  configuration.  In  contrast,  the 
backward-Euler  scheme  is  stable  at  the  Courant  condition.  As  a  result,  the  run  times  using  the 
forward-Euler  scheme  were  about  five  times  as  long  as  those  using  the  backward-Euler  scheme, 
as  the  efficiency  gains  of  the  explicit  forward-Euler  scheme  were  dwarfed  by  the  inefficiency  of 
the  smaller  timesteps.  The  computations  in  the  following  two  sections  were  performed  with  the 
backward-Euler  scheme  only. 


UNCLASSIFIED 


13 


UNCLASSIFIED 


Figure  5.  Plastic  strains  in  a  cylinder  impacting  a  rigid  surface  as  viewed 
from  above.  Part  a  is  the  initial  configuration;  part  b  is  the  response  using 
isotropic  material  constants  in  the  Hill  yield  function;  and  parts  c-h  are  the 
responses  using  50%  increases  in  the  indicated  individual  yield  stresses. 
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3.2  Tensile  Tests 

The  mechanical  response  of  a  tensile  specimen  of  FRP  composite  depends  strongly  on  the 
orientation  of  the  load  with  respect  to  the  principal  material  directions.  When  the  direction  of  the 
load  is  aligned  with  one  of  the  in-plane  principal  directions,  the  fiber  characteristics  dominate  the 
response,  and  the  specimen’s  stress-strain  curve  is  similar  to  that  of  the  fibers.  However,  when 
the  direction  of  the  load  bisects  the  two  in-plane  principal  directions,  the  fibers  allow  much 
larger  strains  before  specimen  failure,  the  matrix  material  contributes  to  the  response,  and  the 
stress-strain  curve  becomes  much  softer  and  more  nonlinear. 

a 


<r 


10  cm 


b 


c  d 


Figure  6.  Tensile-test  results  of  composites  of  S-2  glass  fibers  and  polyester 
matrix  performed  by  Espinosa  et  a!  [5].  The  stress-strain  curves  from  the 
tests  are  shown  in  part  c,  and  from  EPIC  computations  in  part  d. 
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Espinosa  et  al  [5]  perfonned  tensile  tests  on  composites  of  woven  S-2  glass  fibers  in  a 
polyester  matrix.  Part  a  of  Fig.  6  is  a  diagram  of  their  tensile  specimen.  Cut  from  a  panel,  it  is 
10  cm  long,  1.6  cm  wide,  and  0.4  cm  thick.  Loading  is  applied  through  strain-compatible 
fiberglass  end  tabs.  The  angle  9  is  defined  by  the  orientation  of  the  nearest  in-plane  principal 
material  direction  with  respect  to  the  loading  direction.  Specimens  loaded  along  the  fibers  are 
therefore  designated  by  9  —  0°.  EPIC  computations  were  perfonned  for  the  same  values  of  9 
that  were  reported  by  Espinosa  et  al:  9  —  0°,  10°,  20°,  30°,  45°.  Part  b  of  the  figure  shows  the 
mesh  of  3200  hexahedral  elements  that  represents  half  of  the  length  of  the  tensile  specimen.  A 
plane  of  symmetry  was  enforced  at  one  end  of  the  mesh,  and  velocities  were  applied  at  the  other 
end. 

The  stress-strain  curves  reported  by  Espinosa  et  al  are  reproduced  in  Part  c  of  Fig.  6.  With 
increasing  9 ,  the  curves  demonstrate  a  decrease  in  stiffness  and  large  increases  in  both  the 
nonlinearity  of  the  response  and  the  strains  to  failure.  The  material-model  parameters  listed  in 
Table  2  were  used  in  the  EPIC  computations,  with  anisotropy  in  both  the  elasticity  and  plasticity 
models.  These  constants  were  chosen  to  produce  the  stress-strain  curves  shown  in  part  d  of  Fig. 
6.  The  similarity  of  the  stress-strain  curves  in  parts  c  and  d  therefore  does  not  serve  as 
independent  validation  of  the  model  constants,  but  rather  as  a  demonstration  of  the  orthotropic 
model’s  ability  to  simulate  some  of  the  fundamental  characteristics  of  FRP  composites.  These 
characteristics  include  a  nearly  linear  response  to  failure  when  loading  along  a  principal  in-plane 
material  direction,  and  much  greater  plastic  strains  as  the  loading  direction  deviates  from  the 
principal  in-plane  material  direction. 

Table  2.  Orthotropic  Model  Parameters  Used  to  Simulate  the  Tensile  Tests  of  Espinosa  et  al  [5] 


It  should  be  noted  that  the  experimental  results  in  part  c  of  Fig.  6  indicate  failure  of  the  FRP 
composite  at  approximately  the  same  load  when  9  =  0°  and  9  =  10°,  but  at  nearly  twice  the 
strain  when  9  —  10°.  This  behavior  suggests  that  the  failure  is  stress  dependent,  in  contrast  to 
the  plastic-strain-dependent  failure  that  has  been  implemented  in  the  orthotropic  model.  Further 
examination  of  the  failure  characteristics  of  FRP  composites  may  therefore  be  necessary  to  refine 
the  ortho  tropic  failure  model  for  composites. 

Unloading  curves  with  9  —  20°  are  represented  by  dashed  lines  in  parts  c  and  d  of  Fig.  6. 
The  straight-line  path  of  the  computation  reflects  linear  elasticity,  while  the  curved  path  of  the 
experiment  indicates  some  nonlinearity  in  the  elasticity.  Since  the  stress-strain  curves  with 
0  =  0°  are  linear,  it  appears  that  the  strong  experimental  nonlinearities  that  occur  when  9  >  0° 
are  partially  elastic,  and  that  the  nonlinear  elasticity  is  due  to  the  polyester  matrix. 
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3.3  Ballistic  Impacts 

The  computations  in  this  section  demonstrate  the  use  of  the  orthotropic  model  to  evaluate 
FRP  composites  under  ballistic  impact.  The  top  of  Fig.  7  shows  the  plane  of  symmetry  of  the 
impact  configuration.  The  square  target  is  24  cm  wide  and  1.6  cm  thick.  It  is  composed  of  a 
composite  material,  with  the  first  and  second  principal  material  directions  in  the  plane  of  the 
target,  as  indicated  by  the  material  axes  at  the  top  of  the  figure.  The  cylindrical  projectile  is  2.4 
cm  in  diameter  and  4.8  cm  long,  and  it  is  modeled  as  an  isotropic  elastic-plastic  tool  steel  [12]. 


Figure  7.  Plane-of-symmetry  view  (top)  of  the  initial  configuration  of  a 
cylinder  impacting  a  composite  target,  and  the  ballistic  limits  resulting 
from  variations  in  the  individual  yield  strengths  in  the  composite  (bottom). 

Table  3  lists  the  reference  orthotropic  material  parameters  used  in  the  target.  Although  they 
do  not  represent  a  specific  composite  material,  they  include  two  features  of  balanced-weave 
FRPs:  equal  yield  stresses  and  stiffnesses  in  the  in-plane  principal  directions,  and  equal  yield 
stresses  and  stiffnesses  in  the  transverse  shear  directions.  They  also  include  a  transverse  normal 
yield  stress,  a33,  that  is  significantly  less  than  the  in-plane  values. 
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Table  3.  Reference  Material  Parameters  Representing  a  Balanced-Weave  Composite,  and  Used  to 
Compute  a  Ballistic  Limit  of  360  m/s  When  Impacted  by  a  Cylinder  2.4  cm  in  Diameter 


The  computations  were  perfonned  with  an  erosion  strain  of  1.5,  so  that  elements  on  the 
contact  interfaces  are  removed  from  the  computation  when  their  effective  plastic  strains  reach 
150%,  and  the  timestep  does  not  become  prohibitively  small.  Given  the  critical  plastic  strains  in 
Table  3,  the  erosion  strain  of  1.5  ensures  that  elements  will  only  be  discarded  after  they  have  lost 
strength. 

With  the  reference  parameters  in  Table  3,  a  ballistic  limit  of  360  m/s  was  computed.  To 
evaluate  the  influence  of  the  yield  stresses,  ballistic  limits  were  also  computed  using  the 
reference  parameters  with  variations  in  the  individual  yield  stresses.  Variations  of  the  in-plane 
nonnal  yield  stresses,  <JX1  and  o22,  were  coupled,  as  were  variations  in  the  out-of-plane  shear 
yield  stresses,  ay3  and  a23.  For  every  combination  of  material  parameters,  computations  were 
perfonned  at  4  m/s  intervals  to  determine  the  ballistic  limit. 

The  plot  of  ballistic  limits  as  a  function  of  the  yield  stresses  at  the  bottom  of  Fig.  7 
summarizes  the  results  of  the  computations.  The  data  points  from  variations  in  the  same  yield 
stress  are  connected  by  a  line.  One  data  point  in  each  of  these  sets  represents  the  reference 
material  parameters,  and  it  is  indicated  by  a  larger  circle  than  the  other  data  points.  The  four 
large  circles  in  the  figure  therefore  represent  the  ballistic  limit  (360  m/s)  from  the  same 
parameters,  while  the  smaller  circles  each  represent  the  ballistic  limit  from  a  unique  combination 
of  parameters,  with  a  variation  from  the  reference  parameters  in  one  of  the  yield  stresses. 

As  a  yield  stress  increases,  the  ballistic  limit  either  increases  or  remains  constant.  When  the 
ballistic  limit  increases,  the  yield  stress  is  limiting  the  ballistic  performance  of  the  target.  And 
when  the  ballistic  limit  remains  constant,  a  different  yield  stress  is  limiting  the  ballistic 
perfonnance.  In  the  case  of  these  reference  parameters,  it  appears  that  the  nonnal  yield  stresses 
are  limiting  the  ballistic  performance,  but  the  shear  yield  stresses  are  not.  In  addition,  it  appears 
that  the  in-plane  shear  yield  stress,  oy2 ,  will  never  limit  the  ballistic  performance.  Given  the  very 
low  values  of  oy2  for  which  the  ballistic  limit  has  been  computed,  this  observation  may  apply  to 
many  impact  configurations  and  many  FRPs.  It  is  also  consistent  with  the  good  perfonnance  of 
ballistic  fabrics,  for  which  ay2  is  essentially  zero. 

When  the  reference  parameters  are  changed,  or  the  ratio  of  target  thickness  to  projectile 
diameter  is  increased,  the  limiting  yield  stresses  will  likely  change.  In  particular,  the  transverse 
shear  yield  stresses  are  expected  to  limit  a  thicker  target’s  performance,  since  tests  indicate  that 
delamination  becomes  a  critical  mode  of  inelastic  deformation  in  FRP  targets  subjected  to 
ballistic  impacts. 
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As  mentioned  previously,  the  quadratic  form  of  the  Hill  yield  surface  is  not  sufficiently 
general  to  admit  any  combination  of  normal  yield  stresses.  Walker  and  Thacker  [13]  noted  that 
for  orthotropic  plates  the  out-of-plane  nonnal  yield  stress  must  be  greater  than  half  the  in-plane 
nonnal  yield  stress,  2 crj3  >  a\x  —  er^2,  and  proposed  a  quartic  yield  surface  to  circumvent  this 
limitation.  Tensile  tests  of  FRP  composites  generally  register  a  much  smaller  yield  strength  in 
the  transverse  direction,  governed  by  the  matrix  strength  and  delamination,  than  they  do  in  the 
in-plane  principal  directions.  The  Hill  yield  surface  will  therefore  fail  to  capture  the  low 
transverse  tensile  strengths  of  FRP  composites.  The  significance  of  this  shortcoming  in 
computations  of  ballistic  impacts  has  yet  to  be  determined. 


UNCLASSIFIED 


19 


UNCLASSIFIED 


20 


UNCLASSIFIED 


UNCLASSIFIED 


4.0  Summary  and  Conclusions 

This  report  documents  the  development  of  a  material  model  in  EPIC  that  is  orthotropic  in  its 
elasticity,  plasticity  and  failure.  The  model  has  been  developed  to  simulate  the  loading  to  failure 
of  fiber- reinforced  polymer  composites.  Plasticity  is  characterized  by  the  Hill  yield  function, 
and  failure  is  indicated  when  any  component  of  the  material’s  plastic  natural-strain  tensor 
reaches  its  user-supplied  critical  value. 

The  elastic-plastic  stress  update  was  implemented  using  both  forward-Euler  and  backward- 
Euler  finite-difference  schemes  in  time.  Although  the  forward-Euler  scheme  is  explicit,  example 
computations  uncovered  a  maximum  timestep  that  is  significantly  less  than  the  maximum 
timestep  required  by  the  equations  of  motion  -  i.e.,  the  Courant  condition.  On  the  other  hand, 
the  iterative  backward-Euler  scheme  was  stable  in  all  of  the  example  computations  when  the 
timestep  was  governed  by  the  Courant  condition.  In  addition,  an  initial  estimate  of  the  solution 
to  the  iterative  scheme  was  found  that  proved  reliable  and  efficient.  As  a  result,  the  backward- 
Euler  scheme  has  demonstrated  superior  performance,  and  it  has  been  implemented  in  the 
production  version  of  EPIC. 

Computations  were  performed  to  demonstrate  the  effects  of  the  model  parameters,  and  for 
an  initial  assessment  of  the  model’s  applicability  to  composites.  Tensile  tests  on  FRP 
composites  demonstrate  nearly  linear  behavior  to  failure  when  the  direction  of  loading  is  aligned 
with  one  of  the  principal  material  directions,  but  much  greater  strains  to  failure  and  nonlinearity 
when  the  direction  of  loading  deviates  from  the  principal  material  directions.  Computations  with 
the  new  orthotropic  model  were  shown  to  reproduce  this  behavior  when  the  model  constants 
were  chosen  to  reflect  greater  stiffness  and  less  ductility  in  the  principal  material  directions. 

Computations  of  the  ballistic  limits  of  orthotropic  targets  were  included  to  demonstrate  the 
model’s  abilities  and  potential  shortcomings  in  ballistic  applications.  The  computations 
identified  the  components  of  yield  stress  that  limit  the  ballistic  performance  of  the  example 
configuration  and  material  parameters.  An  additional  observation  was  that  the  Hill  yield 
function  requires  that  no  normal  yield  stress  be  less  than  half  the  other  normal  yield  stresses. 
This  restriction  may  pose  a  limitation  to  the  model  for  some  impact  configurations  and 
composite  materials.  However,  in  many  instances  the  transverse  shear  stresses  allow 
delamination  before  the  transverse  normal  yield  stress  is  reached.  As  a  result,  the  practical 
significance  of  the  restriction  in  the  disparity  of  normal  yield  stresses  has  yet  to  be  determined. 
A  second  potential  shortcoming  of  the  Hill  yield  function  is  the  equivalence  of  tensile  and 
compressive  yield  stresses. 

As  more  laboratory  and  ballistic  test  data  become  available,  further  evaluation  of  the  new 
orthotropic  model  will  be  possible.  This  may  lead  to  modifications  or  refinements  in  its  form. 
Potential  improvements  include:  modification  of  the  form  of  the  failure  criteria;  addition  of  an 
optional  yield  surface;  and  transformation  to  a  Total  Lagrangian  formulation,  with  its  precise 
treatment  of  orthotropy  in  the  presence  of  finite  shearing  deformations.  Finally,  the  insight 
gained  from  meso-scale  computations  can  also  be  used  to  improve  the  homogenized  orthotropic 
model,  or  to  substitute  for  missing  test  data  when  orthotropic  model  parameters  are  sought  for 
specific  materials. 
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